Pour cette visualisation, nous avons mobilisé des données provenant du satellite VIIRS (Visible Infrared Imaging Radiometer Suite). Les données téléchargées depuis le service FIRMS de la NASA portent uniquement sur les douze derniers mois. Elles incluent plusieurs variables, parmi lesquelles nous utilisons principalement la variable ‘frp’ (Fire Radiative Power) pour représenter l’intensité des incendies. Cette variable, exprimée en mégawatts (MW), quantifie le taux d’émission radiative des feux actifs, offrant une indication directe de l’intensité du feu. De plus, nous avons défini une zone d’intérêt de 20km autour de l’aire protégée et nous visualisons les données à l’intérieur de cette zone. La carte finale montre ces incendies avec une gradation de couleurs pour refléter l’intensité de chaque incendie.En cliquant sur les points, on peut consulter la date à laquelle le foyer a été enregistré par le satellite, ainsi que la probabilité qu’il s’agisse téellement d’un feu (classé en “élevée”, “nominale” et “faible”).
Code
# Ci-dessous les librairies R utilisées pour cette analyselibrary(tidyverse) # pour la manipulation et visualisation de donnéeslibrary(sf) # pour traiter les données spatialeslibrary(tmap) # pour générer des carteslibrary(aws.s3) # Pour accéder aux données volumineuses stockées en lignelibrary(geodata) # Pour avoir les frontières administrativeslibrary(curl) # Pour téléchargelibrary(lubridate) # Pour manipuler des dateslibrary(httr)library(mapme.biodiversity)library(units)library(arrow)
Code
# Create a dataframe of the coordinatescoordinates <-data.frame(location =c("Beroroha", "Beronono", "Tsivoko", "Makaykely"),latitude =c("21°40.504’S", "21°21.669’S", "21°17.712’S", "21°28.074’S"),longitude =c("45°9.571’E", "45°14.885’E", "45°22.732’E", "45°21.896’E"))# Function to parse DMS and convert to decimal degreesddm_to_decimal <-function(dms){# Extraire les degrés, minutes et secondes parts <-regmatches(dms, gregexpr("[0-9.]+", dms))[[1]] degree <-as.numeric(parts[1]) minute <-as.numeric(parts[2]) /60 direction <-ifelse(grepl("S|W", dms), -1, 1)# Converstion au format décimal decimal <- direction * (degree + minute)return(decimal)}# Apply function to each coordinatecoordinates$latitude_decimal <-sapply(coordinates$latitude, ddm_to_decimal)coordinates$longitude_decimal <-sapply(coordinates$longitude, ddm_to_decimal)# Convertir les coordonnées en objet sfcoordinates_sf <-st_as_sf(coordinates, coords =c("longitude_decimal", "latitude_decimal"), crs =4326)makay_ap <-st_read("data/Makay_wgs84.geojson", quiet =TRUE) %>%rename(Statut = SOURCETHM) %>%mutate(Statut =recode(Statut, "zone tampon"="Zone tampon"))# Fusion des différentes airesmakay_union <- makay_ap %>%st_union() %>%st_sf() %>%st_make_valid()# Define the buffer size (in kilometers) for the square expansionbuffer_size_km <-20# Create a 20km buffer around Makay PAmakay_buffer <- makay_union %>%st_buffer(dist =set_units(buffer_size_km, km)) %>%st_as_sf()# Manual request at https://firms.modaps.eosdis.nasa.gov/download/create.php# Uses the following coordinates : 44.4,-22.4,46.4,-19.9# For VIIRS NOAA-20firms_focus <-read_csv("firms_manual.csv") %>%st_as_sf(coords =c("longitude", "latitude"), crs =4326) %>%filter(lengths(st_within(., makay_buffer)) >0)breaks_frp <-c(min(firms_focus$frp), median(firms_focus$frp), max(firms_focus$frp))breaks_frp_5 <-quantile(firms_focus$frp, probs =seq(0, 1, by =0.25), na.rm =TRUE)tmap_mode("view")tm_shape(firms_focus) +tm_dots(col ="frp", size =1, alpha =0.2,palette =c("blue", "cyan", "yellow", "orange", "red"),border.col =NULL, breaks = breaks_frp_5,popup.vars =c("acq_date", "confidence", "frp"),title ="Fire radiative power") +tm_layout(title ="Feux détectés autour du Malay") +tm_shape(makay_union) +tm_borders(col ="darkblue") +tm_shape(coordinates_sf) +tm_dots(col ="black", size =0.1, labels = coordinates_sf$location) +tm_view(set.view =10) +tm_text(text ="location", ymod =-0.5) +tm_basemap("OpenStreetMap") +tm_scale_bar()
4.2 Feux détectés sur deux décennies
Ces données s’appuient sur les données MODIS qui sont moins précises mais couvrent une période plus longue. On compare le Makay avec une zone périphérique de 60km. This analysis is based on the example produced by Darius Goergen on the Mapme website.
Code
s3_output <-"diffusion/fire-db2.parquet" makay_union <-st_make_valid(makay_union)if (is.na(st_crs(makay_union))) st_crs(makay_union) <-4326# Utilise un CRS métrique pour les buffers (Makay est vers ~45E, 20S -> UTM 38S)makay_m <-st_transform(makay_union, 32738) # EPSG:32738 (UTM 38S)buf20_m <-st_buffer(makay_m, 20000) # 20 km# Donut = buffer - polygonedonut_m <-st_make_valid(st_difference(buf20_m, st_buffer(makay_m, 0)))# Reviens en WGS84 pour les jointures avec (x,y)makay_wgs <-st_transform(makay_m, 4326)donut_wgs <-st_transform(donut_m, 4326)# 1) Ouvrir le dataset Parquet (Arrow détecte les partitions Hive tout seul)ds <-open_dataset(paste0("s3://projet-betsaka/", s3_output), format ="parquet")# 2) Pré-filtrer côté Arrow par bbox + variable (réduction d’E/S)bb <-st_bbox(st_union(makay_wgs))tbl <- ds %>%filter( variable =="fire_occurrences_1year", x >=!!bb["xmin"], x <=!!bb["xmax"], y >=!!bb["ymin"], y <=!!bb["ymax"]# , year >= 2002 # décommente si tu veux mimer l’exclusion 2000–2001 ) %>%select(pixelid, x, y, year, value) %>%collect()# 3) Conversion en points sf puis tag "zone"pts <-st_as_sf(tbl, coords =c("x","y"), crs =4326, remove =FALSE)in_makay <-as.logical(st_intersects(pts, makay_wgs, sparse =FALSE)[,1])in_donut <-as.logical(st_intersects(pts, donut_wgs, sparse =FALSE)[,1])pts$zone <-NA_character_pts$zone[in_makay] <-"Makay"pts$zone[!in_makay & in_donut] <-"Périphérie 20km"pts2 <- pts %>%filter(!is.na(zone))# 4) Agrégats annuels (et total si tu veux)res_by_year <- pts2 %>%st_drop_geometry() %>%group_by(zone, year) %>%summarise(n_cells =n(),sum_fire =sum(value, na.rm =TRUE),mean_fire=mean(value, na.rm =TRUE),.groups ="drop" ) %>%arrange(zone, year)res_total <- pts2 %>%st_drop_geometry() %>%group_by(zone) %>%summarise(n_cells =n(),sum_fire =sum(value, na.rm =TRUE),mean_fire=mean(value, na.rm =TRUE),.groups ="drop" )res_by_year
p_sum <-ggplot(res_by_year, aes(x = year, y = sum_fire, color = zone)) +geom_line() +geom_point(size =1) +labs(x ="Année", y ="Somme annuelle de feux",color ="Zone",title ="Somme annuelle de feux (Makay vs Périphérie 20 km)") +theme_minimal()# 2) Moyenne annuelle par zone (courbe)p_mean <-ggplot(res_by_year, aes(x = year, y = mean_fire, color = zone)) +geom_line() +geom_point(size =1) +labs(x ="Année", y ="Moyenne annuelle de feux par cellule",color ="Zone",title ="Moyenne annuelle de feux par cellule (Makay vs Périphérie 20 km)") +theme_minimal()# 3) Moyenne sur l'ensemble des années (barres)avg_by_zone <- res_by_year %>%group_by(zone) %>%summarise(mean_sum_fire =mean(sum_fire, na.rm =TRUE),mean_mean_fire =mean(mean_fire, na.rm =TRUE),.groups ="drop" )# a) moyenne annuelle de la somme (barres)p_avg_sum <-ggplot(avg_by_zone, aes(x = zone, y = mean_sum_fire, fill = zone)) +geom_col() +labs(x =NULL, y ="Somme annuelle moyenne",title ="Somme annuelle moyenne de feux (2000–2024)") +theme_minimal() +theme(legend.position ="none")# b) moyenne (sur années) de la moyenne par cellule (barres) — optionnelp_avg_mean <-ggplot(avg_by_zone, aes(x = zone, y = mean_mean_fire, fill = zone)) +geom_col() +labs(x =NULL, y ="Moyenne annuelle par cellule",title ="Moyenne annuelle par cellule (2000–2024)") +theme_minimal() +theme(legend.position ="none")# Afficherp_sum
Fire occurence count in Makay PA (source: MODIS)
Code
p_mean
Fire occurence count in Makay PA (source: MODIS)
Code
p_avg_sum
Fire occurence count in Makay PA (source: MODIS)
Les différents ordres de grandeur doivent être interprétés avec précaution car la zone périphérique est plus de deux fois plus vaste que l’aire protégée (3,751 km2 contre 8,249 km2).
